#!/usr/bin/env bash

####################### FUNCTIONS #########################

strand() {
	if [ $1 == A ];then
		all=T
	elif [ $1 == T ];then
		all=A
	elif [ $1 == C ];then
		all=G
	else  
		all=C
	fi
	echo $all
}

############### 

chrpos() {
	echo "chrpos	IDdata	freqdba0	freqdataa0	deltafreq" > "output_liftoversh/informations/keepchrpos_"$name"_"$2".txt"
	echo "chrpos	IDdata	freqdba0	freqdataa0	deltafreq" > "output_liftoversh/informations/outchrpos_"$name"_"$2".txt"
	th=`echo "scale=1;(1/10)"|bc`
	while read line
	do   
		while IFS="	" read chrpos IDdata chr IDdb a0db a1db freqdba0 freqdataa0 deltafreq a0data a1data
			do
				if [ $a0db == $a0data ] && [ $a1db == $a1data ];then
					if test $deltafreq \< $th ;then 
						echo "$chrpos	$IDdata	$freqdba0	$freqdataa0	$deltafreq"  >>"output_liftoversh/informations/keepchrpos_"$name"_"$2".txt"
					else
						echo "$chrpos	$IDdata	$freqdba0	$freqdataa0	$deltafreq" >> "output_liftoversh/informations/outchrpos_"$name"_"$2".txt"
					fi
				elif [ $a0db == $(strand $a0data) ] && [ $a1db == $(strand $a1data) ];then
					if test $deltafreq \< $th ;then 
						echo "$chrpos	$IDdata	$freqdba0	$freqdataa0	$deltafreq"  >>"output_liftoversh/informations/keepchrpos_"$name"_"$2".txt"
					else
						echo "$chrpos	$IDdata	$freqdba0	$freqdataa0	$deltafreq" >> "output_liftoversh/informations/outchrpos_"$name"_"$2".txt"
					fi
				elif [ $a0db == $a1data ] && [ $a1db == $a0data ];then
					newfreqdba0=`echo "scale=2;1-$freqdba0" |bc`
					delta=`echo "scale=2;$newfreqdba0 - $freqdataa0" |bc` 
					newdelta=${delta#-} #absolute value
					if test $newdelta \< $th ;then
						echo "$chrpos	$IDdata	$newfreqdba0	$freqdataa0	$newdelta"  >> "output_liftoversh/informations/keepchrpos_"$name"_"$2".txt"
					else
						echo "$chrpos	$IDdata	$newfreqdba0	$freqdataa0	$newdelta">> "output_liftoversh/informations/outchrpos_"$name"_"$2".txt"
					fi
				elif [ $a0db == $(strand $a1data) ] && [ $a1db == $(strand $a0data) ];then
					newfreqdba0=`echo "scale=2;1-$freqdba0" |bc`
					delta=`echo "scale=2;$newfreqdba0 - $freqdataa0" |bc` 
					newdelta=${delta#-} #absolute value
					if test $newdelta \< $th ;then
						echo "$chrpos	$IDdata	$newfreqdba0	$freqdataa0	$newdelta"  >> "output_liftoversh/informations/keepchrpos_"$name"_"$2".txt"
					else
						echo "$chrpos	$IDdata	$newfreqdba0	$freqdataa0	$newdelta">> "output_liftoversh/informations/outchrpos_"$name"_"$2".txt"
					fi
				else
					echo  "$chrpos	$IDdata	$freqdba0	$freqdataa0	$deltafreq" >> "output_liftoversh/informations/outchrpos_"$name"_"$2".txt"
				fi
			done
	done < $1
	echo -e "see list of SNPs to keep : output_liftoversh/informations/keepchrpos_"$name"_"$2".txt" 
	echo -e "see list of SNPs to del : output_liftoversh/informations/outchrpos_"$name"_"$2".txt" 
}

########################### MAIN #############################

###check for parameters' number
if [ $# -lt 4 ];then
	echo -e "\n ERROR : liftover.sh <date_of_freqfiles> <mergefile> <databasefile> <file1> <file2> ... \n file1 = liftover_wellconverted_dataname1.txt, file2 = liftover_wellconverted_dataname2.txt etc... \n command example : ./liftover4.sh 7212011 ../../db/RsMergeArch.bcp ../../db/NEW_1000/EUR.chrfulldata.db liftover_wellconverted_*\n" 1>&2
	exit 1
fi

### check for frequencies files date
if [ -d $1 ];then # check existence of freqfile
	echo "The first parameter must be a date number(e.g. frequencies_dataname_1232011.out, date_of_freqfiles : 	1232011)"
	exit 2
else
	if [ $(($1+0)) -eq 0 ];then # eq 0 if it isn't a number
		echo "The first parameter must be a date number(e.g. frequencies_dataname_1232011.out, date_of_freqfiles : 	1232011)"
	exit 2
	else
		datefreqfile=$1
	fi
fi
shift

datefile=`date +%m%e%Y` # date mmddyyyy

### mergefile
echo -e "\n************* CHOOSE OR LOAD MERGE FILE *************"
find merge_[0-9]*.mg > /dev/null 2>&1
if [ $? -ne 0 ];then
	echo "load $1"
	cat $1 | tr '	' ' ' | sort -t' ' -k1,1 > "merge_"$datefile".mg"
	mergefile="merge_"$datefile".mg"
else
	echo -e "\nmerge file(s) in your directory : "
	find merge_[0-9]*.mg
	echo -e "\nenter merge filename to use or press \"c\" to load merge parameter's file"
	read choice1
	if [ "$choice1" == "c" ];then
		echo "load $1"
		cat $1 | tr '	' ' ' | sort -t' ' -k1,1 > "merge_"$datefile".mg" 
		mergefile="merge_"$datefile".mg"
	else
		mergefile=$choice1
	fi
fi
shift

### databasefile
echo -e "\n*********** CHOOSE OR LOAD DATABASE FILE ************"
dbname=$1 # save for perl script argument
find db_[0-9]*.db > /dev/null 2>&1
if [ $? -ne 0 ];then
	echo "load $1"
	cat $1 |  awk -F" " '{print $1"-"$3,$1,$2,$4,$5}' | sort -t' ' -k1,1 > "db_"$datefile".db" #sort 1000G database file by chromosome and position
	db="db_"$datefile".db"
else
	echo -e "\ndatabase(s) file(s) in your directory : "
	find db_[0-9]*.db 
	echo -e "\nenter database filename to use or press \"c\" to load database's parameter's file"
	read choice2
	if [ "$choice2" == "c" ];then
		echo "load $1"
		cat $1 |  awk -F" " '{print $1"-"$3,$1,$2,$4,$5}' | sort -t' ' -k1,1 > "db_"$datefile".db" #sort 1000G database file by chromosome and position
		db="db_"$datefile".db"
	else
		db=$choice2
	fi
fi
shift

### create output directories if it doesn't exist yet
if [ ! -d "./output_liftoversh" ]; then
	mkdir output_liftoversh
fi
if [ ! -d "./output_liftoversh/informations" ]; then
	mkdir output_liftoversh/informations
fi

echo -e "\n************** START FILES' TREATMENT ***************\n"
for f in $@;do

	echo "FILE : $f " 
	### datafiles wellconverted "rs" SNPs
	name=`basename $f .txt` # remove ".txt" suffix
	cat $f | sed 's/^chr//' | awk -F"	" '{print $1"-"$2,$4}' | sort -t' ' -k1,1 > "sort_"$name"_"$datefile".txt" # sort by chromosome and position
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB" > "output_liftoversh/informations/summary_"$name"_"$datefile".txt" # header
	join -t' ' -1 1 -2 1 "sort_"$name"_"$datefile".txt" $db | tr ' ' '	' >> "output_liftoversh/informations/summary_"$name"_"$datefile".txt"  # join files by identical chromosome-position
	echo "see info : output_liftoversh/informations/summary_"$name"_"$datefile".txt"

	# select lines with rsID in 1000genomes database 
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB" > "output_liftoversh/informations/rs_"$name"_"$datefile".txt"  # header 
	cat "output_liftoversh/informations/summary_"$name"_"$datefile".txt" | egrep "	[0-9]+	rs" >> "output_liftoversh/informations/rs_"$name"_"$datefile".txt" 

	echo "see info : output_liftoversh/informations/rs_"$name"_"$datefile".txt"

	### datafiles wellconverted "chr-pos" SNPs and frequencies' files
	if [ ! -f "frequencies"${name#liftover_wellconverted}".bim_"$datefreqfile".out" ];then
		echo "for frequencies_dataname.bim_dateoffreqfile.out liftover file must be : liftover_wellconverted_dataname.txt and liftover_failed_dataname.txt and don't forget the first parameter : dateoffreqfile"
	exit 3
	fi
	cat "frequencies"${name#liftover_wellconverted}".bim_"$datefreqfile".out" | sort -t'	' -k2,2 | awk -F"	" '{print $2,$6,$7,$8,$4,$5}' | tr '	' ' ' > tmpfreq.txt
	# select lines with chrposID in 1000genomes database 
	cat "output_liftoversh/informations/summary_"$name"_"$datefile".txt" | egrep -v "	[0-9]+	rs" | sort -t'	' -k2,2 | tr '	' ' ' > "tmpchrpos_"$name"_"$datefile".txt" 
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB	freqDB_a0 freqdata_a0	deltafreq	a0data	a1data" > "output_liftoversh/informations/chrpos_"$name"_"$datefile".txt" #header of chrpos file
	join -1 2  -2 1  "tmpchrpos_"$name"_"$datefile".txt" tmpfreq.txt | awk -F" " '{print $2,$1,$3,$4,$5,$6,$7,$8,$9,$10,$11}' | tr ' ' '	' >> "output_liftoversh/informations/chrpos_"$name"_"$datefile".txt"  

	echo "see info : output_liftoversh/informations/chrpos_"$name"_"$datefile".txt"
	
	cp "output_liftoversh/informations/chrpos_"$name"_"$datefile".txt"  "copychrpos.txt"
	grep "	NA	" "copychrpos.txt" | cut -d"	" -f2 > "output_liftoversh/informations/outNAchrpos_"$name"_"$datefile".txt"
	echo "see list of SNPs to del : output_liftoversh/informations/outNAchrpos_"$name"_"$datefile".txt"
	grep -v "	NA	" "output_liftoversh/informations/chrpos_"$name"_"$datefile".txt" > "modified_chrpos_"$name"_"$datefile".txt" #del lines with "NA" fields
	chrpos "modified_chrpos_"$name"_"$datefile".txt" $datefile #create two files : one of SNPs to del and one of SNPs to keep	

	### datafiles failed
	endfailedfile=`basename ${f#liftover_wellconverted} .txt`
	failedfile="liftover_failed"$endfailedfile
	cat $failedfile".txt" | sed '/#Deleted in new/d' | cut -d"	" -f4 | uniq > "output_liftoversh/informations/exclude_"$failedfile"_"$datefile".txt"
	echo "see list of SNPs to del : output_liftoversh/informations/exclude_"$failedfile"_"$datefile".txt"

	### check RsMerge
	# delete "rs" prefix 
	cat "output_liftoversh/informations/rs_"$name"_"$datefile".txt" | sed ':z;s/rs//;tz' > temporary.txt
	# search rsID 1000G in RsMerge
	cut -d"	" -f4 temporary.txt | sort > listotest.txt
	join -1 1 -2 1 listotest.txt $mergefile > "jointmp2_"$name"_"$datefile".txt"
	cat "jointmp2_"$name"_"$datefile".txt" | sort -t' ' -k2,2 > "jointmp3_"$name"_"$datefile".txt"
	cut -d"	" -f2 temporary.txt | sort > listotest2.txt
	join -1 1 -2 2 listotest2.txt "jointmp3_"$name"_"$datefile".txt" | tr ' ' '	' > "tmpkeep_"$name"_"$datefile".txt"
	echo "rsIDdata	rsIDdb" > "output_liftoversh/informations/keeprs_"$name"_"$datefile".txt"
	cat "tmpkeep_"$name"_"$datefile".txt" | awk '{print "rs"$1,"rs"$2}' >> "output_liftoversh/informations/keeprs_"$name"_"$datefile".txt"
	echo -e "see list of SNPs to keep : output_liftoversh/informations/keeprs_"$name"_"$datefile".txt \n" #-e allows \n interpretation

	###create a list with all SNPs to exclude
	if [ ! -f "exclude"${name#liftover_wellconverted}".bim_"$datefreqfile".out" ];then
		echo "file : exclude_dataname.bim_dateoffreqfile.out doesn't exist"
	exit 3
	else
		cut -d"	" -f2 "exclude"${name#liftover_wellconverted}".bim_"$datefreqfile".out" > "tmpsnpToDel_"$name"_"$datefile".txt"
		cat "output_liftoversh/informations/exclude_"$failedfile"_"$datefile".txt" >> "tmpsnpToDel_"$name"_"$datefile".txt"
		cut -d" " -f1  "output_liftoversh/informations/keeprs_"$name"_"$datefile".txt" > tmpsnpToDel_1.txt
		cut -d"	" -f2 "output_liftoversh/informations/keepchrpos_"$name"_"$datefile".txt" >> tmpsnpToDel_1.txt
		cut -d"	" -f4 $f > tmpsnpToDel_2.txt
		sort  tmpsnpToDel_2.txt tmpsnpToDel_1.txt  | uniq -u > tmpsnpToDel_3.txt
		sort  tmpsnpToDel_2.txt tmpsnpToDel_3.txt | uniq -d >> "tmpsnpToDel_"$name"_"$datefile".txt" 
		cat "tmpsnpToDel_"$name"_"$datefile".txt" | uniq |  grep "^rs" > "output_liftoversh/snpToDel_"$name"_"$datefile".txt"
		echo -e "see summary list of all SNPs to del : output_liftoversh/snpToDel_"$name"_"$datefile".txt \n" #-e allows \n interpretation
	fi

done

echo -e "************* CREATE NEWBIM FILES TO IMPUTE **************"
perl create_bimfile_to_impute.pl $dbname $datefile newbimfile_*$datefreqfile*

echo -e "\n*********************** ALL DONE *************************"


# delete temporary files
rm sort_*.txt 
rm tmpfreq.txt
rm tmpchrpos_*
rm modified*.txt
rm temporary.txt
rm listotest.txt
rm listotest2.txt
rm jointmp*
rm tmpkeep_*
rm tmpsnpToDel_*
rm copychrpos.txt

# command example : ./liftover4.sh 7212011 ../../db/RsMergeArch.bcp ../../db/NEW_1000/EUR.chrfulldata.db liftover_wellconverted_*
